function [D] = Ddavid_chi2_d(X, Y)

M = size(X, 1);
N = size(Y, 1);
D = zeros(M, N);

for i = 1:N
    d = bsxfun(@minus, X, Y(i, :));
    s = bsxfun(@plus, X, Y(i, :));
    D(:, i) = sum(d.^2 ./ (s / 2 + eps), 2);
end
